subroutine srccsd(ie, fij, fab, fia, & 
          iijkl, iijka, iabij, &
          iiabj, iabci, iabcd, &
          t1, t2, r1, r2, &
          fae, fmi, fme, &
          wmnij, wmbej, icore, & 
          maxcore, escf, eccsd) 
use mod_iop
use mod_size
use mod_orbit
use mod_ioff
implicit none
include 'omp_lib.h'
!
integer prec2, isd2, igpu2, totgpu2
common/precsoc/prec2, isd2, igpu2, totgpu2
!
integer, intent(in) :: maxcore 
real*8, intent(inout) :: escf, eccsd
real*4, intent(inout) :: ie(*), fij(*), fab(*), fia(*) 
real*4, intent(inout) :: t1(*), t2(*), r1(*), r2(*) 
real*4, intent(inout) :: iijka(*), iabij(*), iijkl(*)
real*4, intent(inout) :: fae(*), fmi(*), fme(*)
real*4, intent(inout) :: wmnij(*), wmbej(*)
real*4, intent(inout) :: iiabj(*), iabci(*), iabcd(*), icore(*)
!
real*4, allocatable,device :: icoreg(:)
!
integer nsize, i, j, a, b, icrsizg
integer maxvar, maxdim, iter
integer i0, i1, i2, i3, i4, i5, i6, i7, i8, i9
integer i10, i11, i12
real*8 dtmp, errdiis, eccsd2, eccsd3
real*8 eso, tt(20), rt(20), t1diag
real*8 Eccsdio(1000)
!
real*8, external :: ddot, dnrm2, nrm2so
integer, external :: idsymsoc, irpdso
!
zorb(2:8) = 0
if(totgpu==1) then 
   allocate(icoreg(igpusize))
endif
!ldebug = .false.
!ldebug = .true.
!
call myicopyso(ie, ie((nva+noa)*isd+1), nva+noa)
call myicopyso(iabij, r2, nioi3)
if(ldebug) then 
  call printdnrm2('r2 begin ', r2, 3)
endif
call myicopyso(fia, r1, nifvo)
!
maxvar = nifvo + nioi3 
iter = 0 
if(ldebug) then 
  call printdnrm2('iabij', iabij, 3)
endif
!
!--iteration begins
!
do
iter = iter+1
write(6,*)
write(6,'(A, I4)')'=============in iteration:', iter
!
call cpu_time(tt(1)); call walltimeso(rt(1))
call denomt1(t1, r1, ie, 1)
call denomt2(t2, r2, ie, 1, 1, 1)
if(ldebug) then 
 call printdnrm2('r2before doso ', r2, 3)
 call printdnrm2('t2before doso ', t2, 3)
endif
!call myicopyso(t1, t1(ifvo(2)), nifvoaa)
!call iassymso(t2(ioi3(2)), icore, 1, vrta, vrta, popa, popa, 5)
!call icompso(icore, t2(ioi3(1)), vrta, vrta, 0, 1, popa, popa, 1, 0, 1)
!call myicopyso(t2(ioi3(1)), t2(ioi3(3)), nioi31)
!
maxdim = 5
call dosodiis2(iter, t1, maxvar, maxdim, icore, &
              maxcore+maxvar, errdiis)
call effectt2(t2, t1,  icore, 1.d0, no, nv)
if(ldebug) then 
   call printdnrm2('t2after doso ', t2, 3)
endif
!
!call calcene(t1, t2, fia, iabij, eccsd)
!call calcene2(t1, t2, fia, iabij, eccsd2)

 eccsd3 = 0.d0
 call calcene3(t1, fia, nifvoaa, eccsd3, 2.d0) 
 call calcene3(t2(ioi3(2)), iabij(ioi3(2)), nioi32, eccsd3, 2.d0)
 call isymtrso(iabij(ioi3(2)), icore, vrta, vrta, 1, 1, popa, popa, 1, 0, 1) 
 call calcene3(t2(ioi3(2)), icore,          nioi32, eccsd3, -1.d0)
  
eccsdio(iter)=eccsd3
  write(6,'(A,I4, 2F18.12)')'  CCSD and total energy', iter,  eccsd3, escf+eccsd3
!
call effectt2(t2, t1, icore, -1.d0, no, nv)
if(iter>maxiter) exit
!
call cpu_time(tt(2)); call walltimeso(rt(2))
if(ldebug) then 
  write(6,'(A, F10.3, F10.3)')'====time for energy     ', tt(2)-tt(1), rt(2)-rt(1)
endif
!===========prapere for the begin 
call myicopyso(fab, fae, nifvv)
call myicopyso(fij, fmi, nifoo)
call myicopyso(fia, fme, nifvo)
call myicopyso(iijkl, wmnij, nioi1)
call myicopyso(iiabj, wmbej, nioi4)
call myicopyso(fia, r1, nifvo)
call myicopyso(iabij, r2, nioi3)
call cpu_time(tt(3)); call walltimeso(rt(3))
if(ldebug) then 
  write(6,'(A, F10.3, F10.3)')'====time for prepare    ', tt(3)-tt(2), rt(3)-rt(2)
endif
!
!==============n5 term=====================
!
if(ldebug) then 
  call printdnrm2('iiabj', iiabj, 4)
  call printdnrm2('t1-start' , t1, 22)
  call printdnrm2('r1-start' , r1, 22)
endif
call calciiabj(r1, r2,  t1, &              ! rhf finished
          iiabj, icore)
if(ldebug) then 
  call printdnrm2('r1-iiabj' , r1, 22)
endif 
call cpu_time(tt(4)); call walltimeso(rt(4))
if(ldebug) then 
   write(6,'(A, F10.3, F10.3)')'====time for n5 term    ', tt(4)-tt(3), rt(4)-rt(3)
endif 
!=============iabci term==========
!call printdnrm2('iabci' , iabci, 5)
ioiabci=1
if(ioiabci==1) then 
 call calciabci_n5(r1,r2,t1,t2,fae,wmbej,icore)      ! rhf finished 
elseif(ioiabci==0) then 
 call r1ai_t2efmi_iefam(r1, t2, iabci, icore, -0.5d0)! rhf finished
 call r2abij_t1ei_iabej(r2, t1, iabci, icore, 1.d0)  ! rhf finished 
 call wmbej_t1fj_iefbm(wmbej, iabci, &               ! rhf finished
                   t1, fae, icore, -1.d0, 1.d0) 
endif
igpu=0; igpu2=0
call cpu_time(tt(5)); call walltimeso(rt(5))
if(ldebug) then 
  write(6,'(A, F10.3, F10.3)')'====time for wmbej large', tt(5)-tt(4), rt(5)-rt(4)
  call printdnrm2('r1-t2abci' , r1, 22)
  call printdnrm2('r2-n5term' , r2, 3)
  call printdnrm2('fae', fae, 11)
  call printdnrm2('wmbej', wmbej, 4)
endif
!==============================
!-----terms which are independent 
!=============================
!
!==========abci======================
call cpu_time(tt(6)); call walltimeso(rt(6))
if(totgpu==1) then
  i0=1
  i1 = i0 + nifvo*isd
  i2 = i1 + nioi3*isd
  i3 = i2 + nifvo*isd
  i4 = i3 + nioi3*isd
  i5 = i4 + nvvoo*isd*3
  igpuleft = igpusize - i5
  if(i5>igpusize) goto 200
  igpu=1; igpu2=1
  call myicopy_g(r2,    icoreg(i1), nioi3)
  call myicopy_g(t1,    icoreg(i2), nifvo)
  call myicopy_g(t2,    icoreg(i3), nioi3)
  if(iter<2) write(6,*)'calc iabci large in GPU'
  call calciabcig(icoreg(i1), icoreg(i2), &             !rhf finished
                 icoreg(i3), icoreg(i4), iabci, icore)
  call myicopy_c(icoreg(i1), r2,    nioi3)
else
200 if(totgpu==1) then 
     if(iter<2) write(6,'(A, I25, I15, A)')'gpu memory not enough for iabci ', & 
                               i6, i6/1024/1024*4, 'MB'
    endif
   call calciabci(r1, r2, t1, t2, iabci, icore) ! rhf finished
endif
call cpu_time(tt(7)); call walltimeso(rt(7))
 if(ldebug) then  
   write(6,'(A, F10.3, F10.3)')'====time for iabci      ', tt(7)-tt(6), rt(7)-rt(6)
 endif
  igpu=0; igpu2=0
if(ldebug) then 
call printdnrm2('r1-iabci' , r1, 22)
call printdnrm2('r2-iabci', r2, 3)
endif
!
!==========abcd======================
if(totgpu==1) then
  i0=1
  i1 = i0 + nifvo*isd
  i2 = i1 + nioi3*isd
  i3 = i2 + nifvo*isd
  i4 = i3 + nioi3*isd
  i5 = i4 + nvvoo*isd
  if(i5>igpusize) goto 300
  igpuleft = igpusize-i5
  igpu=1; igpu2=1
  if(iter<2) write(6,*)'calc iabcd large in GPU'
  call myicopy_g(r1,    icoreg(i0), nifvo)
  call myicopy_g(r2,    icoreg(i1), nioi3)
  call myicopy_g(t1,    icoreg(i2), nifvo)
  call myicopy_g(t2,    icoreg(i3), nioi3)
!  write(6,*)'finished copy of iabcd'
  call calciabcd(icoreg(i0), icoreg(i1), icoreg(i2), &   !rhf for d finished 
                 icoreg(i3), iabcd, icoreg(i4))
  call myicopy_c(icoreg(i0), r1,    nifvo)
  call myicopy_c(icoreg(i1), r2,    nioi3)
!  call calciabcd(r1, r2, t1, t2, iabcd, icore)
elseif(totgpu==0) then 
300 if(totgpu==1) then 
     if(iter<2) write(6,'(A, I25, I15, A)')'gpu memory not enough for iabcd ', & 
                               i5, i5/1024/1024*4, 'MB'
    endif
  call calciabcd(r1, r2, t1, t2, iabcd, icore)
endif
igpu=0; igpu2=0
call cpu_time(tt(8)); call walltimeso(rt(8))
if(ldebug) then 
  call printdnrm2('r1-iabcd', r1, 22)
  call printdnrm2('r2-iabcd', r2, 3)
  write(6,'(A, F10.3, F10.3)')'====time for iabcd      ', tt(8)-tt(7), rt(8)-rt(7)
endif
call cpu_time(tt(9)); call walltimeso(rt(9))
!call printdnrm2('r1', r1, 22)
!call printdnrm2('r2', r2, 3)
!
!=======    wmbej   ========================================
!
if(ldebug) then 
  call printdnrm2('r1_b_wmbej' , r1, 22)
endif
if(totgpu==1) then 
  i0=1                  !r1
  i1 = i0 + nifvo*isd   !r2
  i2 = i1 + nioi3*isd   !t1
  i3 = i2 + nifvo*isd   !t2 
  i4 = i3 + nioi3*isd   !iijka
  i5 = i4 + nioi2*isd   !iabij
  i6 = i5 + nioi3*isd   !wmbej
  i7 = i6 + nioi4*isd   !icore 
  i8 = i7 + nvvoo*isd*3   
!  goto 400 
  if(i6+nvvoo*isd*4>igpusize) goto 400
    if(ldebug) then 
       write(6,*)'calc wmbej 2 in GPU'
    endif
if(ldebug) then 
    call printdnrm2('r2', r2, 3)
    call printdnrm2('t2', t2, 3)
    call printdnrm2('iabij', iabij, 3)
endif
    igpu=1; igpu2=1
    call myicopy_g(r1,    icoreg(i0), nifvo)
    call myicopy_g(r2,    icoreg(i1), nioi3)
    call myicopy_g(t1,    icoreg(i2), nifvo)
    call myicopy_g(t2,    icoreg(i3), nioi3)
    call myicopy_g(iijka, icoreg(i4), nioi2)
    call myicopy_g(iabij, icoreg(i5), nioi3)
    call myicopy_g(wmbej, icoreg(i6), nioi4)
    call calcwmbej(icoreg(i4), icoreg(i5), &                    ! rhf finished
          icoreg(i2), icoreg(i3), icoreg(i0), icoreg(i1), &
          icoreg(i6), icoreg(i7))
!   call calcwmbejg(icoreg(i4), icoreg(i5), &                    ! rhf finished
!         icoreg(i2), icoreg(i3), icoreg(i0), icoreg(i1), &
!         icoreg(i6), icoreg(i7), wmbej, icore)
    call myicopy_c(icoreg(i0), r1,    nifvo)
    call myicopy_c(icoreg(i1), r2,    nioi3)
else 
400 if(totgpu==1) then 
    if(iter<2) write(6,'(A, I125, I15, A)')'gpu memory not enough for wmbej 2 ', & 
                               i8, i8/1024/1024*4, 'MB'
    endif
!    call printdnrm2('iijka', iijka, 2)
    call calcwmbej(iijka, iabij, &                                ! rhf finished
                t1, t2, r1, r2, &
                wmbej, icore)
endif
  igpu=0; igpu2=0
call cpu_time(tt(10)); call walltimeso(rt(10))
if(ldebug) then 
  call printdnrm2('r1' , r1, 22)
  call printdnrm2('r2-wmbej', r2, 3)
  write(6,'(A, F10.3, F10.3)')'====time for wmbej      ', tt(10)-tt(9), rt(10)-rt(9)
endif
!==============left============================

if(totgpu==1) then 
  i0=1                  !r1
  i1 = i0 + nifvo*isd   !r2
  i2 = i1 + nioi3*isd   !t1
  i3 = i2 + nifvo*isd   !t2 
  i4 = i3 + nioi3*isd   !iijka
  i5 = i4 + nioi2*isd   !iabij
  i6 = i5 + nioi3*isd   !iiabj
  i7 = i6 
  i8 = i7 + nioi1*isd   !fae 
  i9 = i8 + nifvv*isd   !fmi
  i10= i9 + nifoo*isd   !fme
  i11= i10+ nifvo*isd   !icore
  i12= i11+ nvvoo*3*isd
  if(i12>igpusize) goto 500
    if(iter<2) write(6,*)'calc leftterm in GPU'
  igpu=1; igpu2=1
  call myicopy_g(r1,    icoreg(i0), nifvo)
  call myicopy_g(r2,    icoreg(i1), nioi3)
  call myicopy_g(t1,    icoreg(i2), nifvo)
  call myicopy_g(t2,    icoreg(i3), nioi3)
  call myicopy_g(iijka, icoreg(i4), nioi2)
  call myicopy_g(iabij, icoreg(i5), nioi3)
  call myicopy_g(wmnij, icoreg(i7), nioi1)
  call myicopy_g(fae, icoreg(i8), nifvv)
  call myicopy_g(fmi, icoreg(i9), nifoo)
  call myicopy_g(fme, icoreg(i10), nifvo)
  call calcsmall(icoreg(i4), icoreg(i5), &
            icoreg(i2), icoreg(i3), icoreg(i0), icoreg(i1), &
            icoreg(i8), icoreg(i9), icoreg(i10), &
            icoreg(i7), icoreg(i11))
  call myicopy_c(icoreg(i0), r1,    nifvo)
  call myicopy_c(icoreg(i1), r2,    nioi3)

else
500 if(totgpu==1) then 
     if(iter<2) write(6,'(A, I125, I15, A)')'gpu memory not enough for left term ', &
                               i8, i8/1024/1024*4, 'MB'
    endif
 
call calcsmall(iijka, iabij, &            ! rhf finished
          t1, t2, r1, r2, &
          fae, fmi, fme, &
          wmnij, icore)
endif
igpu=0; igpu2=0
call cpu_time(tt(11)); call walltimeso(rt(11))
if(ldebug) then 
   write(6,'(A, F10.3, F10.3)')'====time for left       ', tt(11)-tt(10), rt(11)-rt(10)
   call printdnrm2('r1 end', r1, 22)
   call printdnrm2('r2 end', r2, 3)
endif
!
write(6,'(A, F10.3, F10.3)')'====time for iteration  ',  rt(11)-rt(1)
!
enddo
write(6,*)'the final energy is', Eccsdio(iter)
open(555,file='CCSD_ENERGY')
write(555,*) ESCF, Eccsdio(iter)
close(555)

      WRITE(6,99)
99    FORMAT(T3,'     Summary of iterative solution of CC equations ')
      WRITE(6,100)
100   FORMAT(6X,59('-'))
      WRITE(6,101)
101   FORMAT(T24,'Correlation',T46,'Total',/,T8,'Iteration',T26, &
             'Energy',T46,'Energy')
      WRITE(6,100)
      DO I=1,ITER
       WRITE(6,1000) I-1, Eccsdio(I), Eccsdio(I)+ESCF, ' DIIS '
1000   FORMAT(T10,I4,T19,F18.12,T39,F19.12,T59,A6)
      ENDDO
      WRITE(6,100)

      T1DIAG = 0.D0
      T1DIAG = nrm2so(irpdso(1,vrta, popa, 1), T1, 1)**2
      T1DIAG = DSQRT(T1DIAG/dble(NOA))
      WRITE(6,*) '     T1 Diagnostic =', T1DIAG
!
      WRITE(6,103) ESCF
      WRITE(6,102) ECCSD3 + ESCF
103   FORMAT('        E(SCF)            = ',F20.12)
102   FORMAT('        E(CCSD)           = ',F20.12)


!
!do i = 1, iter
!   write(6,'(A, I4, F18.12)')'ECCSD with ITER', i, Eccsdio(i)
!enddo
!
return
end 
